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NOMENCLATURE 


Symbol 


What the symbol represents 


Savin [ 1 ] 


u 


jk 


-Q- 


R, R-, 


R, 


eq. (2.1.1) 

E 

G 

£. 


n 

r 

X 

r 

r . 


U 


ci‘ 


N 

1 

g. h 
G, H 
U, P 
I 
0 


Name of the author and the reference number 

Stress tensor 

Displacement component 

Traction component 

Boundary 

Domain 

The variable above the bar appears is prescribed 

as a boundary condition 

Error terms on the domain and different parts of 

the boundary 

Weighting field 
Equation number 
Strain tensor 
Dira Delta function 
Young’ s Modulus 
Rigidity Modulus 
Poisson’s ratio 
Kronecker Delta 

Outward normal to the boundary at x 
Source point 
Field point 

Distance between x and 
Directional distance between x and 
Derivative of r with respect to x. 

Coefficient of ui when i is a point on the 

boundary 

Number of boundary elements 

Length of the element ^ ^ 

Submatrices obtained after integrating u^ . & p^ . 
Assembled matrix of g and h 

Column vector of displacement and traction 
’Unit’ column vector 
’Zero’ column vector 

Stress in the direction B (hoop stress) 

Applied stress 


continued 
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Symbol 


What the symbol represents 


A, B, C 
d 

a , b 
w 


s ll' s 12? 

22' s 66 } 


The three points at which stress concentration 

factors have been computed 
Diameter of the circular cutout 
Minor and major axis of the elliptical cutout 
Width of the finite strip and finite plate 
Young’s modulus in longitudinal and transverse 
direction for composite materials 
Compliance matrix components for 
Composite materials 


Consistent tensor notation has been used throughout. 
Subscripts i,j,k and 1 have been used as indices, the range 
being 1 to 2, unless otherwise indicated. Repetition of the 
index indicates the summation of the variable over the range 
of the index and a ’ , ’ represents differentiation with 
respect to that index of the variable. 



ABSTRACT 

Boundary element method has been used to compute the 
stress concentration factors for various plane stress 
problems, with circular and elliptical cutouts. The problems 
include an infinitely long strip, but finite in its width and 
a finite plate, both of them under uniform tensile streeses . 
The media selected are both isotropic and orthotropic in 
nature. Most of the problems have been done using both 
linear and constant elements. Wherever possible, results 
have been compared with experimental values and analytical 
solutions . 
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CHAPTER 1 


INTRODUCTION AND LITERATURE SURVEY 


loJL 


The boundary element method (BEM) , appropriately called 
the boundary integral equation method (BIEM) is being 
increasingly used to solve solid mechanics problems 
especially when it concerns the stresses on the boundary or a 
crack tip, for example, stress concentration factor, fracture 
mechanics and contact problems. In this work, an effort has 
been made to study the stress concentration factors for 
various classical and other related problems using BEM. 

For elastostatic problems, the method typically starts 
with weighted residual statement of the equilibrium equation. 
The weighting functions are chosen in such a way that the 
error distribution terms over the domain become zero and the 
"domain terms disappear"; only the boundary integrals will 
remain. With discretization of the boundary and substitution 
of proper boundary conditions, namely, displacements and 
tractions, this boundary integral equation can be solved to 
obtain the unknown values on the boundary, which can be later 
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utilised to compute the domain properties ( stresses, 
displacements etc.) 

la. 2 Literature Sur.Y.ez. 

The problem of an infinite plate with a circular cutout 
has been solved analytically by Kirsch, which has beeen 
referred to by most of the literature in solid mechanics, 
including Savin[lJ. On the same lines many problems have 
been solved and results have been presented by Savin[l]. The 
problems include a finite strip (infinite in its length and 
finite in its width) with symmetrically placed circular hole. 
Heywood[2] has done photoelastic analysis of the same 
problem and his results match with that presented by 
Savin [ 1] . 

Later on finite element method was used by Chong and 
Pinter [3] to solve the same problem; their results are 
slightly different from that predicted by Heywood[2] , as d/w 
(d is the diameter of the circular hole and w is the width of 
the infinitely long strip) approaches 1. But for all the 
other values of d/w, their results compound the validity of 
the results presented by others [1,23. 
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Based on the guidelines provided by Brebbia [4,5,6], 
Boundary Element Method was used to solve the same problem; 
though, with constant element the results were not agreeing 
very well with others [1,2,3], linear element gave the results 
which were very close to that presented by them[l,2,3]. 
Using the program developed for constant element, various 
problems of finite strip with elliptical hole and finite 
plate with circular and elliptical hole were solved. 

In the orthotropic case, Lekhnitskii[7] has solved a 
few problems of elliptical and circular hole in an infinite 
medium, as a plane stress problem. Rizzo and Shippy[8] , 
later presented fundamental solutions for plane stress 
orthotropic problems and a procedure exactly same as that of 
isotropic boundary element was used by them to solve a few 
problems in this medium. Their results were agreeing very 
well with that presented by Lekhnitskii [7 ] . Mahajerin and 
Sikarskie[9] have used almost a similar scheme to predict 
the stress concentration factor in a double lap joint under 
tension . 

The problems of elliptical hole in an infinite 
orthotropic medium were solved, based on the fundamental 
solution presented by Rizzo and Shippy[8] , otherwise the 
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method is very similar to that used for isotropic medium. 
The results obtained were matching within 10% of those given 
in Lekhnitskii[7] . This program of constant element was used 
to get the data regading the stress concentration factors for 
finite strip and finite plate with circular and elliptical 
cutout . 

Results have been presented for finite and infinite 
plate with circular and elliptical hole for both isotropic 
and orthotropic media. Most of the cases were studied with 
both linear and constant boundary elements. Wherever 
possible results have been compared with analytical and 
experimental results. 



CHAPTER 2 


BEM FORMULATION 


2uA 


Statement of the problem in elastostatics without body 
force is as follows [ 4 ] : 

Ot, . = 0 with boundary conditions 

JK, J 


u k 

ii 

on 

n 


p k 

ii 

»d 

on 

r 2 

(2.1.1) 


where is the stress tensor 

u? is the displacement component 
p£ is the traction component 

_ means the variable is prescribed as a boundary 

condition 

I - is the boundary enclosing the domain _Q_ 

Tensor notation has been used throughout; repetition of 
an index denotes summation over the range of that index and a 
’ , ’ indicates partial differentiation with respect to the 
varimbl«^ n( j ex (here, for example X j ) . The numbers within the 
square brackets indicate the reference number. 

The generalised weighted residual statement used in 
finite elements and other techniques is [4] 
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( 2 . 1 . 2 ) 


| R w dil 


r R 0 w dr - Jr. (aw/an) dr 


The same weighted residual statement written for the 
elastostatics problem defined by eq. (2.1.1) will be[5] 


/°jk, j u k dj:r= / (p k _p k )u k dr ■ J <u k' a k )p k dr 


where indicates the weighting field. 


Integrating eq. (2.1.3) once 


/ 0 jw € jk d - a - = - I p k u k dr “ I p k u k dr + / <a k' u k )p ^ dr 

z - 'a. m Cj 

where k indicates the strain tensor. 

Noting that [4] Oj k £* k = cr\ k € jk , we have 

j^k £ jk da= -/ p k u k dr -f p k a k dr + =[ <n k- u k )p k dr 


Integrating once again 


(2.1.3) 


(2.2.1) 


/°jk,j u k da = -/Vk^ -/ p k u k dr + l n k p J dr + ,J u k p k dr 

\x M ‘ { '2. 

The above equation without boundary condition can be 


written as 


/°jk, j u k dJi = JVk dr -/ p k u k dr 


( 2 . 2 . 2 ) 
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Now we try to eliminate the domain integral by 
satisfying the equilibrium equation in the following 
manner [6] 

°jk,3 + = 0 (2 - 2 - 3) 

where A * is the Dirac Delta function and hence represents a 
unit force at a point i in direction 1. A traction and 
displacement field denoted as and can be obtained 
from this equation; these will be the fundamental solutions 
to be used in eq. (2.2.2). Substituting eq. (2.2.3) into 
eq. (2.2.2), we obtain 


+ i u k p !k dr = s 


p k u ik dr 


(2.2.4) 


% 

where u^ and p^ are fundamental displacement and traction 


in direction k due to unit force in direction 1. 



Equation (2.2.2) has been solved analytically and the 
fundamental solutions have been presented in many of the 
text-books [4 , 5 , 6] . The fundamental solution to be used in 
the present formulation are that of two-dimensional isotropic 
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plane strain problem, which are as follows: 


u 


* 

i j 



[(3-4V)ln (-)£,, + r r -I 

87TG( 1 -V) r 1J ’ 1 > i-' 

- 1 — <?r 

— {(l-2*)£, . + 2r r } 

47T( 1-V) r — c3n 1J ’ A ,J 


-(1- 2V)(r 


. n .-r 
i J 



(2.3.1) 


Here, the point at which the unit force is applied is 
called the source point (denoted as p ) and the point at which 
we compute the numerical value of the fundamental solution is 
called the field point (denoted as x) . The field point is 
always located on the boundary, whereas the source point may 
be inside the domain or on the boundary. 


In the above fundamental solution and elsewhere 

G is the rigidity modulus 
v is the Poisson’s ratio 

r is the distance between the source point and 

the field point 

n . is the outward normal to the boundary at the 
J field point 

8 . . is the Kronecker delta 
U 

All the derivatives are taken with respect to the field 
point x. 
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To obtain fundamental solution for the plane stress 
case , the substitution of by >* /( 1+ >> ) will suffice . 


2LA Boundary Integral Equation 


In eq. (2.2.4) the point i is inside the domain. We 
can take it to the boundary as a limiting case and this 
results in a slightly different equation of the form[5] 



(2.4.1) 


This is the boundary integral equation (BIE) and can be 
solved numerically to obtain the unknown boundary values (u 
on and p on Tg ) • Once all the boundary values are 
known, eq. (2.2.4) can be used to obtain displacement ( s) in 
the domain. Alternately, we may wish to obtain stresses in 
the domain directly by using 


2G .Ju, 

s . ± + g 


1J 1-2? 1J £x, 


au. 

o)x . 
J 


r$U 




where (o>u./o>x.) can be obtained from 
making the substitution for COu^/Ox^.) 


(2.4.2) 

eq. (2.2.4). After 
eq. (2.4.2) becomes 


CT . 
1J 



P k di- 


Is 

r 


ki3 U k dP 


where 


(2.4.3) 
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4a( 1-y) r 


(l-2V){S kir>j+ S kJ r i i-£ ij r_ k } + 2r > .r >J r |kj 


s, . . 

ki j 


2ff{ 1-V) r 


2-{a-2V)£ u r k ^(S lk r + £. k r a ) 
— d>n 


“ 4r ,i r , j r ,k> +2 ''' (n i r ,j r ,k +n j r ,i r ,k : 
.(l-2W(2n k r ]i r J +n.£. k «.S jt ) 




(2.4.4) 
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CHAPTER 3 


NUMERICAL IMPLEMENTATION 


Er Qc edure 


Writing eq. (2.4.1) for any point on the boundary, 

have : 


c . .u . 
1J 3 


I p . u . 
J 1 3 3 


.dr 


r 



(3.1.1) 


The main steps involved in numerical implementation can 
be written down as: 


(i) Discretization of the boundary {“ into N 
elements over which the displacements and the 
tractions may be assumed to be piecewise 
interpolated. This way eq. (3.1.1) gives 2 linear 
algebraic equations with 4N vairables. 

(ii) Similar equations for all i = l,N 
will be 2N in number, which means it will give us 
a system of linear algebraic equations with 2N 
unknowns. For each node the number of degrees of 
freedom is 4, namely, Ug, P^ an <* Pg, out of 


we 
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which 2 will be specified as boundary conditions 
and the rest of the 2 unknowns per node will make 
it to 2N unknowns . 

(iii) Use any standard method to solve this system 
of equations and obtain the other boundary values. 

(iv) The same procedure is adopted while computing 
the stresses or displacements in the domain, except 
that these values will be obtained directly after 
computing their respective boundary integrals. 

3+2 IzBsa ol Element 

The boundary elements, for two-dimensional problems, as 
the name indicates, will be line elements on the boundary. 
Usually the degrees of freedom will be traction (p) and 
displacement (u) in x and y directions, numbering to 4. 
Depending upon the distribution of these properties over the 
element, the elements can be classified as constant, linear, 
quadratic and so on. 

For constant element, the distribution of u and p will 
be uniform or constant and the representative point will be 
the geometrical mid-point of the element itself. For linear 
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element, a linear distribution is assumed between the two 
extreme ends of the element and for quadratic element a 
parabolic variation may be assumed with another point in the 
middle of the element (refer fig. 3.2.1). 

In this work, most of the computations have been done 
using constant element and, unless otherwise indicated, all 
the numerical implementations have been done for the 
constant element. 


3llA PlagrstlaatiQn and Computation ol En*. 


Equation (3.1.1) in a discretized form can be rewritten 
as follows: 


c . .u . + 
1 3 3 




(3.3.1) 


Integration over each element j is done by using 
Gaussian Quadrature. The number of points taken for this is 
usually 4; in some cases variable number of points have been 
used. The (directional) distances (r.^) and other derivatives 
are given by the following expressions (see fig. 3.3.1) 
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FIG. 3.3.1 


Integration on the boundary 




FIG. 3.3.2 

Analytical evaluation of 


FIG. 3.7.1 
ds/ dx^ 



r i = V x > - x i‘p> 

r = ^ r i r i jl/2 
r . = (Sr/3x. (x) ) = r./r 

i JL JL X 

(Sr /an) = r .n. 

,:L 1 (3.3.2) 

As we see the numerical computation of u. . and p. . 

10 

over j with respect to i does not pose any problem when 

The values obtained after numerically computing the integral 

of u. , and p. . are called as g and h submatrices, each of 
ij U 

which is 2x2. But when i=j, u. . has ln( 1/r) singularity 
which can be evaluated analytically [5] (refer fig. 3. 3. 2), the 
results of which have been reproduced below: 


g ll ~ 


1 

87tG(l-J3) 


1 

(3-4V)(l-ln-) + 
2 



g 12 " g 21 


r l r 2 


8*G(1-V)1 


g 


22 


8rcG( 


[73- 

1-)P) u 


4>3) (l-ln~) + - 

2 \ 1 


f] 


(3.3.3) 


where 1 is the length of the element, 
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Hence diagonal g submatrices can be evaluated 
analytically. Even though the diagonal h submatrices have not 
been computed we will start assembling them in accordance 
with eq. (3.3.1). 

LA Assembly ol ills EiAsr si lasd Eaaa.il.Qii U L LA 1 

For any point i on the boundary, eq. (3.3.1), after 
discretization and assembly will take the following form: 


: ii u i + [b 


h il b i2 


‘i 5 • 


l iJ 


u ] 

u r 


u 


u 


N 


[« 


fi il s i2 


s 


u *• 




P J 


■N 


Varying i from 1 to N we have the following system 
of equations, represented in the matrix form. 


17 



hi hz 

j 21 .22 

;; .2 j 

{jlN 

h 2N 


U 1 

^2 

**il hz 

:: N 

.... h iN 



_^N1 h N2 

•• h Nj 

h NN_^ 


__ U N_ 





N 



fll 

g 12 ' ' 

g l j 

g 1N 


Pi 

f21 

g 22 ' • 

?2j 

g 2N 


Po 

: 2 

fil 

fi2 • ' 

fij 

g iN 

4 

t * 

__5ni 

g N2 • ■ 

g Nj 

g NN 


L P N_ 


where h^ = 

We have evaluated neither nor c^ 

notation eq. (3.4.1) can be rewritten as 

[H] {U} = [G] {P} 


(3.4.1) 


In matrix 

(3.4.2) 


But now h^ can be directly obtained from the consideration 
of rigid-body movements [4] in a bounded or unbounded body. 
If we have a unit displacement in any one direction eq. 
(3.4.2) then becomes 

[H] {I 1 } = {0} 

in which 1^ is a vector defining unit rigid displacement in 
direction 1. Therefore the diagonal submatrices are simply 
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( 3 . 4 . 3 ) 


h 


ii 


S h 


ij 


Thus there 
explicitly, 
(infinite or 


ii 


was no necessity to find out 
But the same condition for an 
semi-infinite) becomes[4] 

= I 


Hi 1J 


°ii and b ii 
unbounded body 

(3.4.4) 


SalYing ills system 


Now eq. (3.4.2) is complete in the sense that we know H 
and G matrix completely, and we can directly substitute the 
boundary conditions in the form of displacements and/or 
tractions, at each node. The total number of unknowns after 
this will be 2N and eq. (3.4.2) can be rearranged as: 

[A] {X} = {B} (3.5.1) 

where A contains the rearranged coefficients, B is the 
column-vector (obtained by multiplying the proper 
coefficients with the known boundary conditions) and X is 
the unknown column-vector. X in eq. (3.5.1) can be obtained 
using standard procedures available for solving system of 
linear algebraic equations. 
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To keep track of the known/unknown tractions and 
displacements, a system of coding has been used. For 

example, K0DE=1 means traction specified and K0DE=0 means 
displacement prescribed at that node, in x or y direction. 
This coding system will be able to sort out the X vector (now 
containing the solution) properly into traction and 
displacement vectors. 

Domain Displacements and Stresses 

Once we know the displacements and tractions over the 
entire boundary, then, for displacement we can use eq. 

(2.2.4) and for stress eq. (2.4.4) can be used. The 

discretization to be adopted will be the same as that used 
for eq. (3.3.1) and exactly the same procedure will be 

employed for this. But now the only difference will be that, 
there will be no singular integrals, because the point under 
consideration is inside the domain. 

2-JL Surface St rsaasa 

Evaluation of stresses on the boundary may also be done 
taking the load point to the boundary, as a limiting case of 
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equation eq. (2.4.4). But this will lead to very inaccurate 
results, mainly because the tensor components involve 
singular integrals. Hence another scheme presented by Rizzo 
and Shippy[8] can be very easily implemented without having 
to undergo the cumbersome procedure of evaluating singular 
integrals. The scheme is outlined below: 


We have seven unknowns on the boundary, namely 3 stress 

components ( <J ~ . .) and 4 displacement gradient components 

(u. .). Also we have 

1 1 0 


= Q~. ,n. 

1 3 i 


(3.7.1) 


along with stress-strain relationships (3 in number). This 
amounts to 5 equations. We need two more equations. They 
can be got by evaluating du^/ds (along the boundary) 
numerically and equating it to analytical expression for the 
same . 


du. du. dx 1 du. dx 9 

1 X JL ^ X Cm 

ds dx^ ds dXg ds 

= u^ 1 (dx 1 /ds) + u i g ( dxg/ds ) 

From f ig . 3.7.1 


( dx^/ds ) = 

- sin 6 

II 

1 

ro 

(dXg/ds) = 

cos e 

n 


(3.7.2) 


( 3 . 7 . 3 ) 
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Substituting this into eq . (3.7.2) 


(du./ds) = u. 9 n. - u. . n 9 

± - , £ j . 1,1^ ( 37 . 4 ) 

where s has been used to indicate the boundary. 

Now we have a system of 7x7 linear algebraic equations, 
which can be solved to give the stresses on the surface or 
the boundary . 


Co mposition ol 

As we have seen already h^ can be computed by 
considering the rigid-body movements, either in a bounded or 
unbounded medium. But for constant elements, for which the 
boundary is smooth[4], an analytical evaluation of h ii can be 
done. By letting the point i go to the boundry, the first 
boundary integral in eq. (2.4.1) can be evaluated as[5] 



5 0 
-0.5 




Hence submatrix will be 


r i 

0 

' JL i 

- 0.5 0 


0.5 0 

Lo 

1__ 

.■T ' 

; -i 

_0 - 0 . 5 _ 


_0 0 . 5 __ 
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CHAPTER 4 


RESULTS AND DISCUSSION 


4lU. Program 

A simple program with constant boundary element was 
developed to start with, to handle two dimensional 
elastostatic problems in isotropic media with single 
continuous boundary. Then additional features like surface 
stresses, variable Gauss points for integration (depending 
upon the distance from the boundary element and its length) 
were introduced. The program was also extended to handle two 
boundary problems. In fact, it can easily be extended to take 
into consideration any number of boundaries. Using this 
program stress concentration factors were obtained for finite 
and infinite plate with circular and elliptical holes at its 
centre. The same problems were repeated with linear element 
program available in [4] . 

Then a constant element program was developed for 
orthotropic media, with almost the same features as for the 
corresponding isotropic case. The major changes required 
were: fundamental solutions u^ and p^ and their 
derivatives. These have been discussed in detail in 
succeeding sections. This program was used to compute the 
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stress concentration factors in orthotropic media for similar 
problems. Before this a check was made to test the 
authenticity of the program by comparing the results of the 
classical problems presented by Lekhnitskii [7] . In all the 
cases &$/p h as heen computed where p is the applied stress. 

l.i.1 Isotropic Problems 

The classical problem of infinite plate with a hole was 
solved; the resulting stress concentration factor was 3.001 
as against 3 obtained analytically by Kirsch. For the case 
with elliptical hole, a result of 4.962 was obtained, the 
analytical result being 5. In these problems, where constant 
element was used, half of the plate was considered for 
analysis instead of the usual one quarter. The reason for 
doing this was that the constant element’s representative 
point is the mid-point of the element and at points where we 
have corners, we cannot expect to compute the stresses 
accurately (eg. , the point of maximum stress in the problems 
under consideration, that is point B in fig ,4. 2.1); however 
small the element is chosen at this point, the stresses 
computed will not be at point B, but slightly away from it 
and hence it is bound to be inaccurate. A half plate can 
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consider such a point B by having a constant element 

symmetrical about y axis and in this case the mid-point of 

the element itself happens to be B. 

For a d/w ratio of 0.1, the number of elements used were 

54, 19 of them on the semi -circle, The problem was repeated 

for various d/w ratio from 0.2 to 0.9, for a finite strip and 

also for a finite plate with circular and elliptical hole. 

As this ratio is increased, the hole, either circular or 

elliptical also increases in size. Hence, a different and 

very fine mesh size on the boundary was used each time. 

T-fa=ffirgfo fhe mesh size has not been reported for each of the 

case' the validity of the same was checked by observing the 

composition of h-.r submatrix, which has to be [0.5 0 

IP 0.5. 

in the case of constant element, as already indicated 

All the problems were solved using Linear Element 
Program given in [1]. Because of the nature of the linear 
element, the advantage was taken of the bi-axial symmetry, 
and only a quarter of the plate was considered for analysis. 

The results of finite strip with (i) circular hole and 
(ii) elliptical hole under uniform tension at infinity have 
been presented in Table 4.1 and 4.2, respectively. The 
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results of finite plate under uniform tension with ( i ) 
circular hole and (ii) elliptical hole have been presented in 
Table 4.3 and 4.4, respectively. 

For the case of finite strip under uniform tension at 
inifinity, with circular hole, the FEM solution by Chong and 
Pinter[3], the experimental results by Heywood[2] and the 
analytical solution by SavinEl] have been presented for ready 
reference . 

For the maximum stress concentration factor at B, it can 
be seen that the constant element results match with other 
results for lower ratios of d/w and deviate very much as d/w 
approaches 0.9. On the other hand, the linear element 
results almost match even at a higher ratio of 0.8. For the 
point C, where the stress is compressive, the results 
obtained by constant element are very much less than those 
by linear element and are not accurate. This is to be 
expected since the constant element is represented by a point 
which will be just above C. The linear element, however, 
gives a more accurate result, since the nodes are at the 
extreme ends for this type of element, (refer fig. 4.2.2) 
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Stresses obtained for the point A first decrease from 1 
to somewhere near 0.75 and then increases upto 4 as d/w ratio 
increases. Before discussing this, let us have a look at the 
modelling done for this case (refer fig. 4. 2. 3). The model 
almost looks similar to that of a beam loaded in the 
longitudinal direction, with a notch included at its base. 
The boundary conditions of this model are, of course, 
slightly different from that of a beam. From this figure, we 
can say that apart from a tensile stress, there will also be 
a moment which will alter the state of stress at A and B. 
For the point B the bending moment will yield more tensile 
stresses, whereas for the point A it will be compressive. 
Because of this, the overall stress at A tends to reduce upto 
a certain value of d/w. Beyond this d/w ratio, the plate 
does not behave anymore like this and its short area at the 
base, entirely takes up the applied load; with the result the 
stress distribution is non-linear across the base, with both 
the points A and B under larger tensile stresses. 

But from Tables 4.3 and 4.4, which show the results for 
finite plate with circular and elliptical hole, the behaviour 
is very much as we expected. The stresses at A decreases 
from almost p to near 0 for d/w = 0.4, and then becomes 
negative (or compressive). In the case of finite plate, it 
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is obvious that, the model looks like a L-beam for circular 
hole (refer fig. 4. 2. 4). This clearly shows why the stress at 
A goes on decreasing as d/w ratio increases from 0.1 to 0 . 9 . 
But for the case with elliptical hole, once again the effect 
is not that much as that with a circular hole. This may be 
because the ellipse approximates to a notch at the base; but 
the plate is not long enough to behave like a beam. Hence 
the state of sti*ess at A for this case is somewhere between 
that of the previous two cases (Tables 4.2 and 4.3). 

Qrthgtropic Case 

The major changes required for the orthotropic case were 
(i) the fundamental solution, (ii) the stress-strain 
relationship and (iii) the derivatives of the fundamental 
solution which were later used in computing the displacement 
and stresses in the domain. 

The fundamental solution in the case of orthotropic case 
is more complicated than that of isotropic case and involves 
a lot of material dependent constants to be evaluated. The 
plane stress fundamental solution given by Rizzo and 
Shippy[8] has been reproduced here: 
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- ^[^Aglnrj - lnr 2 ] 

U 12 = U 21 = ~ K a A l A 2 (6 l~^2 ) 
u 22 = -^CA^ln^/j^ - A^lnr 2 /v^] 


P U = K a tA l /C ^ r 2 ) ‘ V^l^Vk 

p ?2 = K a [M lV r ? - M 2 A l /r l ] 

4l = K aWl A l /( ”4 r ?) - M 2 A 2 /(^r|)3 

p 22 = W“i r ?> - V'^l^Vk 


where 

= 2 (t V*2 )s 22 

K'l + ^g = ( 2s i2 +s 66 )/s 22 

= s ll /s 22 

A i = s 12 " ^i S 22 

2 2^,2,,. 
r i = X 1 (X 2 /W i ) 

= arctanCxg/Cx^^O ] 
~ jc^ngJST; x 2 n l y/ ^i 


(4.3.1) 


(4.3.2) 


( 4 . 3 . 3 ) 
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The derivatives of the fundamental solutions in x and y 
directions are presented in Appendix A. This along with the 
plane stress constitutive equation was used for solving 
problems in orthotropic media. For all the orthotropic 
problems, the material properties chosen were that of birch 
plywood referenced in Lekhnitskii [7 ] . They are as follows: 

S 11 = 8.333 x 10 ^ sq-cm/kgf 

= -0.6 x 10 ® sq-cm/kgf 

S 22 = -667 x 10 6 sq-cm/kgf 

Sgg = 143 x 10 ^ sq-cm/kgf 

First the problem of an infinite plate with an 
elliptical hole was considered for which analytical solutions 
are available in Lekhnitskii [7] . The results of both 
analytical and BEM solution are given in Table 4.5. 

All the problems were done considering them as two- 
boundary problems, because of which the number of elements 
increased. 36 elements were used on the ellipse inside and 
44 elements on the outer boundary, which was 100 times away 
from the ellipse as was the size of the ellipse. By St. 
Venant’s Principle, the effect of the elliptical hole was not 
considered while imposing the boundary conditions on the 
symmetrical elements of the outer boundary (see fig. 4.3.1 ). 
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On the elliptical hole p = 0 and p =0 were the boundary 

x y 

conditions . 

Table 4.5 clearly shows the accuracy of boundary 
element method for these type of problems in orthotropic 
media as well. Then a series of problems involving finite 
strip and finite plate with circular and elliptical holes 
were solved using the constant element program, the results 
of which have been presented in Tables 4.6 thru 4.11 

Table 4.6 gives the stress concentration factors of a 
finite strip with circular hole, at points A, B and C (see 
the figure above the table) using constant boundary element 
and the medium being mildly orthotropic (for Birch-plywood 
E^/Eg = 2). The direction of loading is along the 

Longitudinal axis of the fibres. The pattern in which the 
stresses at the three mentioned points vary with the sise of 
the hole is very similar to that of the isotropic case shown 
in table 4.1; only the values of the stress concentration 
factors are slightly higher than the corresponding values for 
the isotropic medium. 

For the case of finite strip with elliptical hole, once 
again the same kind of similarity is exhibited as that of 
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isotropic case. For example in isotropic case the stress 
concentration factor at B is more for elliptical hole than 
that of circular hole. The same pattern is seen in 
orthotropic case also 

For finite plate with circular and elliptical hole 
(Tables 4.7 and 4.8) another set of data was generated with 
transverse fibres. As has been reported by Lekhnitskii[7] 
the stress concentration factor (for maximum value of stress) 
for transverse fibres is less than that for longitudinal 
fibres. This feature is seen in both tables 4.7 and 4.8. 
Another behaviour to be noted is that at A stresses become 
compressive for both longitudinal and transverse fibres, like 
that of isotropic medium. 

Another set of data was obtained for galss-epoxy 
composite material and this has been presented in tables 4.9 
to 4.11. The properties of this composite material are [9] 

7.1053 x 10~ 7 sq-cm/kgf 

_7 

-1.4921 x 10 sq-cm/kgf 
67.6423 x 10~ 7 sq-cm/kgf 
167.6848 x 10~ 7 sq-cm/kgf 


S 11 = 
3 12 " 
s 22 
S 66 = 
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The results in tables 4.9 to 4.11 for the case of 
glass-epoxy can be compared with those in tables 4.6 to 4.8 
for birch-plywood. It is easily seen that the maximum stress 
concentration factor (at B) is always greater for glass- 
epoxy, whereas the minimum stress concentration (at C) is 
always less. This is the case with longitudinal fibres only; 
the situation is just the reverse for transverse fibres. 
This pattern in the stress concentration may be attributed to 
the strong orthotropic nature of glass-epoxy, for which 
E^/Eg is 9 as compared with birch-plywood, for which the same 
ratio is 2. 

The stress pattern at A, however, is different 
that at B and C. In the case of glass-epoxy, the stress 
is never compressive, with transverse fibres, while 
situation was just the opposite for birch-plywood. 


from 
at A 
the 
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TABLE 4.1 


E - 2 000 000 kgf/sq-cm 

Point B 

*5/t> 


(BEM) 

d/w Const. Linear Savin Heywood FEM 


0.1 3.001 3.05 3.03 3.032 

0.2 3.039 3.218 3.14 3.14 

0.3 3.239 3.435 3.36 3.347 3.449 

0.4 3.562 3.797 3.74 3.693 3.83 

°- 5 4.07 4.379 4.32 4.25 4.413 

0.7 6.017 6.847 6.747 7.166 

0.8 8.086 9.524 10.002 10.276 

0.9 12.81 15.734 20.01 19.544 



d/w A C (-) 




Const. 

Linear 

Const. 

Linear 

0. 

1 

0. 

.977 

0. 

98 

0. 

, 658 

0. 

. 999 

0. 

.2 

0. 

957 

0. 

949 

0. 

621 

1. 

. 148 

0. 

3 

0. 

909 

0. 

891 

0. 

707 

1. 

,282 

0. 

4 

0. 

891 

0. 

823 

0. 

928 

1. 

.451 

0. 

5 

0. 

892 

0. 

758 

1. 

203 

1. 

.644 

0. 

7 

1. 

24 

0. 

755 

1. 

463 

2. 

.078 

0. 

8 

2. 

25 

1. 

056 

1. 

.725 

2. 

.224 

0. 

9 

6. 

467 

4. 

033 

1. 

959 

2. 

.474 
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TABLE 4.2 


E = 

2 000 

000 kgf/sq- 

cm b/a 

= 2 


» = 0.1 

____ 



°S/f> 




b/w 


A 


B 

Linear 



— 

Const 

Linear 

Const. 

Const. 

C (-) 

Linear 

0 . 1 

0. 985 

0. 99 

4.962 

4.935 

0. 523 

1 . 

0. 2 

0. 977 

0.978 

' 4.913 

5.029 

0. 571 

1.042 

CO 

o 

0. 945 

0. 949 

5.152 

5.263 

0. 620 

1.117 

0.4 

0. 935 

0. 905 

5.537 

5.642 

0.764 

1 . 223 

0. 5 

0. 926 

0. 844 

6.143 

6.245 

0. 941 

1.358 

0.7 

1.008 

0.71 

3.448 

8.836 

1.207 

1.716 

00 

o 

2.061 

0.836 

11.082 

11.921 

1.636 

1.937 

0.9 

6.059 

2.845 

17.02 

19.38 

2.202 

2.36 




TABLE 4.3 


E = 

21 000 

000 kgf/sq- 

•cm 

v~,/f 



V = 0.3 

d/w 

Const . 

A 

Linear 

Const. 

B 

Linear 

Const . 

C (-) 
Linear 

i — 1 

o 

0. 938 

0. 965 

2.88 

3.071 

0.472 

1.057 

0.2 

0. 838 

0.846 

3.106 

3.34 

0.725 

1.272 

0.3 

0. 643 

0.601 

3.533 

3.841 

1.081 

1.69 

0.4 

0. 303 

0. 133 

4.222 

4.705 

1.618 

2.44 

0.5 

-0. 252 

-0.765 

5.313 

6.212 

2.44 

3.789 

0. 6 

-1.076 

-2.58 

7.015 

9.01 

3.649 

6.36 

0.7 

-2.005 

-6.622 

9.513 

14.734 

5.163 

11.7 

0.8 

-1.833 

-16.62 

12.616 

27.882 

7.061 

23.88 

o 

1 

CD 

4.091 

-35.54 

16.48 

54. 

7.23 

+43.8 
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TABLE 4.4 

E = 21 000 000 kgf/sq-cm b/a = 2 


°s/p 



Const . 

21 

Linear 

Const . 

B 

Linear 

Const . 

C (-) 
Linear 

0 . 1 

0.952 

0.978 

4.68 

0.481 

1.036 

1.036 

0.2 

0.891 

0. 905 

4 . 988 

5. 268 

0.703 

1 . 156 

0.3 

0.777 

0.765 

5.461 

5.785 

0.934 

1 . 377 

0.4 

0.589 

0. 525 

6.15 

6.556 

1.238 

1.737 

0.5 

0.303 

0.131 

7.097 

7.66 

1.654 

2.3 

CD 

o 

-0.098 

-0.49 

8.436 

9.257 

2.22 

3.125 

0. 7 

-0. 571 

-1.423 

10.428 

11.672 

2.966 

4.632 

0.8 

-0.798 

-2.65 

13.795 

15.685 

3.895 

6.18 

0.9 

1.442 

-2.98 

21,517 

24.13 

4.928 

8.77 


A 
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TABLE 4.5 


BIRCH PLYWOOD 


b/a = 2 

CONSTANT ELEMENT 



r e/f> 



Major Axis 
inclination 

MAXIMUM 

MINIMUM 

(-) 

to fibres 
( degrees ) 

BEM 

LEKH. 

BEM 

LEKH. 


LONGITUDINAL FIBRES 


90 

10.62 

9.91 

0.669 

0.71 

60 

7.14 

8.2 

0.669 

1.11 

45 

6.099 

5.84 

0.696 

0.46 

30 

5.062 

4.92 

0.677 

0.77 

0 

3.24 

3.23 

0.615 

0.71 


TRANSVERSE FIBRES 



180 

7.56 

7.2 

1.33 

1.41 

150 

5.6 

6.14 

1.25 

1 . 53 

135 

4.67 

4.98 

1.39 

1 . 57 

120 

3.82 

3.8 

1.28 

1.54 

90 

2.56 

2.57 

1.13 

1.41 
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TABLE 4.6 


BIRCH PLYWOOD 


CONSTANT ELEMENT 


°e/p 

d/w or b/w A B C (-) 


CIRCULAR HOLE 


0. 1 

0.996 

5.677 

0.671 

0.2 

0. 971 

5.809 

0.745 

0. 3 

0.942 

6.055 

0.849 

0.4 

0. 925 

6.452 

0.976 

0.5 

0. 937 

7.07 

1.122 

0.6 

1.002 

8.044 

1.287 

0.7 

1.187 

9.664 

1.468 

0.8 

1 . 778 

12.684 

1.654 

0. 9 

4.877 

19.99 

1.863 


ELLIPTICAL 

HOLE 



( b/a = 2 

) 


0.1 

1.003 

10.71 

0.673 

0.2 

0.997 

10.898 

0.71 

0.3 

0.988 

11.247 

0.77 

0.4 

0. 983 

11.804 

0. 854 

0.5 

0.991 

12.666 

0.96 

0.6 

1.034 

14.021 

1.093 

0.7 

1.184 

16.273 

1.26 

0.8 

1.682 

20.493 

1.466 

0.9 

4.312 

30.826 

1.708 


Q 


T J3 

T^N 


w b 
1 


U' 


u-a-»i 


Fiber direction 


43 



TABLE 4.7 


BIRCH PLYWOOD constant element 


d/w 

L 

A 

T 

B 

L 

T 

C (- 
L 

-) 

T 

0. 1 

0.963 

0. 968 

5.922 

4.33 

0.736 

1.443 

0.2 

0.833 

0.871 

6.659 

4.706 

0. 989 

1.873 

0.3 

0. 563 

. 0.692 

1.716 

5.305 

1.468 

2.653 

0.4 

0.082 

0. 412 

9.22 

6.169 

2.254 

3.942 

0.5 

-0. 737 

-0.2 

11 . 366 

7.429 

3.494 

5.988 

CD 

O 

-2.13 

-0.69 

14.72 

9.336 

5.44 

9.308 

0.7 

-4.386 

-1.633 

19.864 

12.26 

8.394 

14.495 

o 

CO 

-7.134 

-2.208 

27.29 

16.427 

12.253 

21.771 

0.9 

-5.347 

1.744 

36.876 

22.633 

15.706 

28.609 


L - Longitudinal Fibres (with usual properties [6] ) 



44 



TABLE 4.8 


BIRCH 

PLYWOOD 


b/a = 2 

V~e/P 


CONSTANT ELEMENT 


A 


B 


C (-) 


b/w 

L 

T 

L 

T 

L 

T 

0. 1 

0.974 

0.98 

11.11 

7 . 762 

0.719 

1 .412 

0.2 

0.887 

0.92 

12.3 

8.293 

0.88 

1.648 

0.3 

0.704 

0. 806 

13.9 

9.078 

1.165 

2.065 

0.4 

0.379 

0.628 

15.893 

10.094 

1.594 

2.684 

0.5 

-0.136 

0. 385 

18.32 

11.376 

2.196 

3.538 

0.6 

-0.879 

0.102 

21.457 

13.045 

3,018 

4.679 

0.7 

-1.809 

-0.132 

25.492 

15.395 

4. 113 

6.159 

0. 8 

-2.59 

-0.024 

31.307 

19.379 

5.504 

7.982 

0.9 

-1.353 

1.985 

43.339 

29.714 

7.11 

10. 169 


A 



Fiber direction 
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TABLE 4.9 


GLASS EPOXY 


°e/> 


CONSTANT ELEMENT 


d/w or b/w 

A 

B 

C (-) 


CIRCULAR 

HOLE 


0. 1 

0.974 

7.146 

0.339 

0.2 

0. 923 

7.23 

0.361 

0.3 

0.912 

7.32 

0.38 

0.4 

0. 899 

7.617 

0.407 

0. 5 

0.884 

8.227 

0.45 

0. 6 

0.889 

9.275 

0. 506 

0.7 

0. 962 

11.137 

0. 566 

0. 8 

1.306 

14.62 

0. 643 

0.9 

3.671 

23.183 

0.723 


ELLIPTICAL 

. HOLE 



(b/a = 

2) 


0.1 

0. 982 

13.666 

0. 333 

0.2 

0.954 

13.774 

0.345 

0.3 

0. 964 

13.885 

0.354 

0.4 

0.971 

14.314 

0.372 

0.5 

0.971 

15.211 

0.403 

0. 6 

0.978 

16.736 

0.449 

0.7 

1.035 

19.426 

0. 505 

0. 8 

1 . 330 

24.289 

0.58 

0.9 

3.379 

36.598 

0.759 



46 



GLASS 

-EPOXY 


TABLE 4. 

°i/f> 

10 

CONSTANT ELEMENT 

d/w 

A 

L 

T 

B 

L 

T 

C 

L 

(-) 

T 

0.1 

0.962 

0.982 

7.18 

2.78 

0.342 

3.031 

0. 2 

0.818 

0.939 

8.328 

2.935 

0.465 

3.728 

0.3 

0.489 

0.887 

9.882 

3.212 

0.695 

4.855 

0.4 

-0.16 

0.831 

12.103 

3.638 

1.062 

6.562 

0.5 

-1.332 

0.794 

15.121 

4.269 

.1.617 

9.077 

0. 6 

-3.37 

0.815 

19.89 

. 5.198 

2.423 

13.183 

0.7 

-6.714 

1.078 

26.878 

6.559 

3.507 

19,438 

0.8 

-11.202 

2. 348 

37.217 

8.56 

4.631 

28.476 

0.9 

-11.744 

7.785 

50.902 

12.811 

5.298 

36.136 



GLASS -EPOXY 


TABLE 4.11 

b/a = 2 


CONSTANT ELEMENT 




b/w 

L 

A 

T 

B 

L 

T 

I 

I o 

1 

1 ^ 

1 

I 

1 

j 

(-) 

T 

0. 1 

0.974 

0.991 

13.726 

4.53 

0.336 

3.034 

0.2 

0.864 

0.971 

15.696 

4.698 

0.425 

3.38 

CO 

o 

0. 62 

0.938 

18.234 

4.99 

0.579 

3.94 

0.4 

0. 145 

0.901 

21.559 

5.433 

0.805 

4.703 

0.5 

-0.696 

0.877 

25.592 

6.084 

1.118 

5.69 

o 

CD 

-2.086 

0.907 

31.191 

7.067 

1.545 

6.976 

0.7 

-4. 19 

1.09 

38.484 

8.668 

2.125 

8.623 

o 

CD 

-6.748 

1.751. 

48.529 

11.733 

2.87 

10.586 

O 

CO 

-6.968 

4.596 

63. 195 

20.238 

3.686 

13.217 


A 



Fiber direction 
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APPENDIX - A 


The derivatives of the fundamental solution in the case 
of plane stress orthotropic problems is given below. These 
are used in computing the stresses in the domain. 


u ll, 1 

= Vjt^A^/r 2 ) 

- (J^A|/r“)] 

U ll, 2 

= V2 [A ? /(r 2^> 

- A 2 /(r 2 ^)] 

U 12 , 1 

= U 21 , 1 = K a A l A 2 x 

2 Cl/(J^r 2 ) - l/(^ d r 2 )] 

u 12 , 2 

= U 21 , 2 ~ K a A l A 2 x 

1 tl/(J^r 2 ) - l/(^r 2 )] 

u 22 , 1 

= K a x 1 [A 2 /(r 2 J^ 1 ) 

- A 2 /(r|y^)3 

U 22 , 2 

= w4 /ir l°<? > 

- 4^44^ 1 

p n,i 

= K a [2x 1 x k n k {A 1 /(^ i r 2 ) - Ag/U^r*)} 


- n 1 {A 1 

/(J^r 2 ) - A 2 /(^r 2 )>] 


PU.Z 1 - A 2 /<«?S» 

- n 2 {A 1 /(^r^) _ A2 / (J ^ r 2 )}] 

p 12 1 = K a ^ A 2 {2x l M l /r l ' J ^ n 2 /r l } " A l {2x l M 2 /r 2 ~ \ n 2 /r l }} 
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P 12 , 2 = K a CA 2 {2x 2 M l / K r l ) + V ( ^i r ? )} 


- A 1 {2x 2 M 2 /U a r^) + ^/(^r 2 )}] 

P 21 , 1 = K a CAiC2x 1 M 1 /rJ - n 2 /r 2 }/ 

- A 2 {2x 1 M 2 /r 2 - n 2 /r 2 >/ ] 

P 21,2 = K a tA 1 { 2 x 2 M 1 /(^r^) + n^C r 2 )}/^ 

- A 2 {2x 2 M 2 /(^r 2 ) + n^Cv^r 2 )}/^] 

p 22 , 1 ’ K a [ 2 x l x k n k {A l /( ' £ k r l ) ” A 2 /( ^ r 2 5 * 


p 22 , 2 


- n^Aj/C^r 2 ) - A 2 /(J^r 2 )}] 
= K a [2x 2 x k n k {A 1 /(°<f i r^) - A 2 /(^' r 2 ) } 

- n^A^G^r 2 ) - A 2 /(^r 2 )}] 


where the usual notations have been used and all the 
constants have been explained in the fundamental solution 
itself . 
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